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Abstract. Jeans showed analytically that, in an infinite uniform-density isothermal gas, plane-wave perturbations collapse to 
dense sheets if their wavelength. A, satisfies A > = {na- jGp^ (where a is the isothermal sound speed and is the 

unperturbed density); in contrast, perturbations with smaller A oscillate about the uniform density state. Here we show that 
Smoothed Particle Hydrodynamics reproduces these results well, even when the diameters of the SPH particles are twice the 
wavelength of the perturbation. Our simulations are performed in 3-D with initially settled (i.e. non-crystalline) distributions 
of particles. Therefore there exists the seed noise for artificial fragmentation, but it does not occur. We conclude that, although 
there may be - as with any numerical scheme - ‘skeletons in the SPH cupboard’, a propensity to fragment artificially is evidently 
not one of them. 
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1. Introduction 

Stars form through the collapse and fragmentation of molec¬ 
ular clouds. This is a highly chaotic and non-linear process; it 
involves many different physical effects, of which the dominant 
one is arguably gravitational fragmentation; and it involves a 
very large dynamic range of physical scales and complex ge¬ 
ometries. Because the process is chaotic and non-linear, nu¬ 
merical simulations have a central role to play in understand¬ 
ing the interplay between the different physical effects. Because 
gravitational fragmentation is a critical effect, it is essential 
that numerical schemes are able to capture this effect prop¬ 
erly, i.e. that true gravitational fragmentation is not suppressed 
by inadequate resolution, and that artificial fragmentation does 
not occur. Because the process involves a large dynamic range 
of physical scales and complex geometries. Smoothed Particle 
Hydrodynamics has been used extensively to simulate star for¬ 
mation (e.g. in 2004 alone, Bonnell, Vine & Bate, 2004; Clark 
& Bonnell, 2004; Delgado-Donate, Clarke & Bate, 2004a; 
Delgado-Donate et ah, 2004b; Goodwin, Whitworth & Ward- 
Thompson, 2004a,b,c; Hennebelle et ah, 2004; Hosking & 
Whitworth, 2004a,b; Jappsen & Klessen, 2004; Li, Mac Low & 
Klessen, 2004; Rice et ah, 2004; Kurosawa et ah, 2004; Price 
& Monaghan, 2004a,b; Schmeja & Klessen, 2004; Whitehouse 
& Bate, 2004). 

In this paper we present a new demonstration of the ability 
of Smoothed Particle Hydrodynamics to simulate gravitational 
fragmentation properly, even at very poor resolution, using the 
plane-wave analysis initially performed by Jeans (1929). In 
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Section 12 we briefly describe SPH in general, and the imple¬ 
mentation we use in particular. In Section|2 we define - from 
first principles - the resolution required for simulating grav¬ 
itational fragmentation, i.e. the so-called Jeans Condition. In 
Section|2 we describe the Jeans Test, and in Section|3 we ex¬ 
plain how the initial conditions for the test are set up. In Section 
12 we present the results of the test, emphasising how poor the 
resolution can be made before the results are significantly cor¬ 
rupted. In Section 0 we derive and demonstrate a correction 
term for use when the Jeans condition is not satisfied or only 
weakly so. In Section|2 we summarise our main conclusions. 

2. Smoothed Particle Hydrodynamics 

Smoothed Particle Hydrodynamics (hereafter SPH) is a 
Lagrangian, particle-based scheme, first formulated by Lucy 
(1977) and by Gingold & Monaghan (1977). The fluid is repre¬ 
sented by an ensemble of particles having masses m,, positions 
r., velocities v,., internal energies m, and smoothing lengths h.. 
Apart from its gravity, the influence of particle i extends only 
to radius r = 2h., and is weighted by the smoothing kernel 
Wirjhff There need be no grids or symmetry constraints, and 
mass is conserved automatically, so there is no need to solve a 
continuity equation. Local functions of state can be evaluated 
at an arbitrary position r by summing contributions from all 
the particles j whose smoothing kernel overlaps r, weighted by 
W{\r-r \lh ). During the evolution, v,. is updated using a sum of 
terms representing hydrostatic, viscous, gravitational and mag¬ 
netic accelerations. Similarly m, can be updated using a sum of 
terms representing adiabatic, viscous and radiative ‘heating’. In 
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principle, the radiative heating can be coupled to a treatment of 
radiation transport, but in practise this is normally very compu¬ 
tationally expensive, and often a barotropic equation of state is 
used instead. 

The SPH code we use here is dragon. This is an exten¬ 
sively tested code (Goodwin, Whitworth & Ward-Thompson, 
2003a), which uses an octal spatial decomposition tree (Barnes 
& Hut, 1986) to speed up the calculation of gravitational ac¬ 
celerations and to find lists of neighbours. Gravity is kernel- 
softened. Particle smoothing lengths, are adapted to give 
A/fjEiB = 50 neighbours. A 2nd-order Runge-Kutta time- 
integration scheme is used, with multiple particle time-steps 
regulated by a Courant-like condition. The code uses time- 
dependent viscosity with = 0.1 (Morris & Monaghan 
1997); the option exists to reduce the effective shear viscos¬ 
ity still further using the Balsara switch (Balsara 1995), but 
this option is not exercised here. Periodic boundary conditions 
can be imposed. If, as here, self-gravity is involved, the Ewald 
method is used (Hernquist et ak, 1991; Klessen, 1997). Since 
we are imposing an isothermal equation of state (i.e. P - a^p) 
the energy equation need not be solved. 

3. Resolution requirements for gravitational 
fragmentation 

The resolution requirements for numerical simulations of grav¬ 
itational fragmentation were first addressed systematically by 
Truelove et al. (1997) and Bate (& Burkert (1977). Truelove et 
al. (1997) used an Adaptive Mesh Refinement (AMR) Finite 
Difference (FD) code to show that grid-based codes must use 
cell sizes d satisfying a Jeans Condition of the form d < 
/1jeans/ 4, where Tjeans - is the local Jeans length, 

a is the local isothermal sound speed, and p is the local den¬ 
sity. If this Jeans Condition is not met in FD simulations, artifi¬ 
cial fragmentation can occur, particularly if artificial viscosity 
is used. 

In this context, the configuration initially explored by Boss 
& Bodenheimer (1979) has acquired the status of a standard 
test, since it is believed that the simulations of this config¬ 
uration presented by Truelove et al. (1998) achieved conver¬ 
gence. The initial configuration involves a spherical cloud with 
mass M - IM^, radius R = 5 x 10'® cm, uniform density 
p - 3.8 X 10^'^gcm^^, uniform temperature T = lOK and 
uniform angular speed Q. = 7.2 x 10 (hence ratio of ther¬ 
mal to gravitational energy a - 0.26 and ratio of rotational to 
gravitational energy (3 = 0.16). An azimuthal m - 2 density 
perturbation with fractional amplitude A = 0.1 is then imposed 
and the subsequent isothermal evolution is followed. In the 
Truelove et al. (1998) AMR simulation of this configuration, 
a binary system forms, with a filament between the two com¬ 
ponents, and - as predicted by Inutsuka & Miyama (1992) - 
the filament does not fragment, but rather collapses to a singu¬ 
larity. This is in direct contrast with the earlier FD simulations 
of the Boss & Bodenheimer configuration reported by Burkert 
& Bodenheimer (1993). They used nested grids to increase the 
central resolution, but did not satisfy the Jeans Condition, and 
this resulted in fragmentation of the filament (into regularly 
spaced condensations with planetary masses). 


The corresponding Jeans Condition for SPH has been de¬ 
rived by Bate & Burkert (1997), who show that SPH models 
gravitational fragmentation properly (i.e. artificial fragmenta¬ 
tion is avoided and true fragmentation captured) provided (i) 
the gravity softening and particle smoothing have similar scales 
(this is achieved automatically here with kernel-softened grav¬ 
ity), and (ii) the local Jeans mass is resolved at all times. They 
suggest that the minimum mass that can be resolved by SPH 
is given by where Af^Em ~ 50 is the mean 

number of neighbours and m is the mass of a single SPH par¬ 
ticle (here assumed to be universal). In a subsequent paper 
(Bate, Bonnell & Bromm, 2002), this has been revised down 
to = I.SA/^eib^’ therefore we have elected to put 

Mmin = AffjEiB ^ to use the factor by which Mje^ns exceeds 
as a measure of the accuracy of the code. Requiring 


M., 


= hi ffi < A4 
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or equivalently 

^JEANS > 


( 3 ) 


i.e. the Jeans length should exceed the diameter of an SPH par¬ 
ticle. 

Whitworth (1998) has shown analytically that with a stan¬ 
dard smoothing kernel, and a gravitational softening length 
comparable to the kernel smoothing length, the only gravita¬ 
tional condensations which can form in SPH must (a) be gen¬ 
uinely gravitational unstable, and (b) be - at least approxi¬ 
mately - resolved. Thus, failing to satisfy the Jeans Condition 
simply suppresses true fragmentation, rather than promoting 
artificial fragmentation. The present paper confirms these re¬ 
sults. 

Kitsionas & Whitworth (2002) have shown that standard 
SPH simulates the Boss & Bodenheimer test as well as AMR, 
provided sufficient particles are used to satisfy the Jeans 
Condition. Moreover, the number of particles can be greatly 
reduced (and with it the amount of memory and computation 
required) by implementing on-the-fly particle splitting. Despite 
the transient high spatial-frequency noise introduced by parti¬ 
cle splitting, the filament between the two binary components 
shows no sign of fragmenting, even when the simulation is 
evolved to densities almost 100 times higher than Truelove et 
al. (1998) achieved, and the computational cost is greatly re¬ 
duced. Thus, ~ 600,000 particles are required for standard 
SPH to follow this test with an isothermal equation of state 
to p = 6 X 10“*gcm“^ (a more stringent test than attempted 
by Truelove et al., 1998). With particle splitting this can be 
achieved with ~ 45,000 particles initially, and fewer than 
~ 145,000 at the end, with huge savings in the required CPU. 


4. The Jeans Test 

Although the Boss & Bodenheimer test is a demanding one, 
it does not have an analytic solution, established beyond all 
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reasonable doubt, and therefore it is appropriate to explore a 
test which does have an analytic solution. 

Consider a static infinite medium, with uniform density p^, 
and uniform and constant isothermal sound speed a, and as¬ 
sume that, even though it is self-gravitating, it is in equilibrium. 
This assumption is usually referred to as the ‘Jeans swindle’, 
since the medium cannot strictly be in self-gravitating equilib¬ 
rium if the density is finite. However, in the simulations pre¬ 
sented below the uniform density unperturbed state is effec¬ 
tively in equilibrium, in the sense that it can be evolved for a 
long time without changing. 

Now suppose that the medium is perturbed so that —> 
p„ -I- Pj and v„ = 0 — » v„ H- Vj = Vj. The continuity, Euler and 
Poisson equations reduce to the forms 


<9Pi 

dt 

= -Po^-V, 

(4) 

5v, 

Vp 

(5) 

dt 

=-^ 

Po 

VY, 

= 47rGp| 

(6) 


where <p^ is the gravitational potential due to the perturbed den¬ 
sity (e.g. Binney & Tremaine, 1987). Eliminating Vj and 
from Eqns. 0 to (|^ then yields 



Fig. 1. The solid line represents one wavelength of the imposed 
sinusoidal plane-wave perturbation. The other lines represent 
the smoothing kernels used in the results displayed in Eig. |2j 
7? = 0.5 (dotted line) the very well resolved case; 7? = 1.0 
(dash-dot line) the marginally resolved case; and 7? = 2.0 
(dashed line) the under-resolved case. The kernels are all scaled 
so that the integrated area under the kernel is equal to the area 
under the perturbation. 


^ -47rGp„p. = 0. (7) 

Substituting a plane wave of the form p^ (r, t) - A gi(kx±o>t) 
Eqn. m gives the dispersion relation 

Lo^ - - AnGpg. ( 8 ) 


Hence we can identify a critical Jeans wave-number, 

, (47rGp„)'/2 


(9) 


and correspondingly a critical Jeans wavelength, 


2 TT 


71 a 
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,2\1/2 


( 10 ) 


is real and the perturbation oscillates. Taking the real parts of 
Eqns. (II 1> and \\2\ . 


P,(r,f) 

v,(r,f) 


- A p^ cos 


A oj, 

-i sin 

k 


l2nx\ 


^ 1 cos(w,f) , 

(14) 

(2nx\ 


1 ^ lsin(w,/)e^. 

(15) 


and the oscillation period is 



n 

Gpo 


1/2 


(T2 - T2)E2 

'• IPAWC 


(16) 


Eor long wavelength perturbations (T > k < 

the dispersion relation (Eqn. 0 indicates that is negative, 
and therefore a>^ is imaginary. Defining 


separating short wavelength perturbations (which oscillate) 
from long wavelength perturbations (which grow). 

To set up an initially stationary plane-wave perturbation, 
we superimpose two plane waves of equal amplitude and wave¬ 
length, travelling in opposite directions: 


7, 


2n a 


1 



■JEANS 


^y/2 

d2j 
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we can put Then, taking the real parts of Eqns. (E) 

and o, we have 
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i 

^ 2n x'\ 
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yi(kx-oip) ^i(kx+ojp)'^ 

(11) 

p,(r,f) = Ap„ cos I 

K 4 J 

^^iikx-iop) _ ^iikx+ojp)'^ - 

(12) 

v.(r,f) = ^ sin| 

|2 7r.r^ 

1 sinh(y,f) 


(18) 

(19) 


Eor short wavelength perturbations (A < k > 
the dispersion relation (Eqn.0 indicates that (Y is positive, and 
therefore (switching from k to A), 


CO, 




1/2 


(13) 


The time for the perturbed density on the plane x = 0 to grow 
from Ap„ to cosh(l) Ap„ = 1.54Ap„ is 

1 

47rGp„ 



1/2 


(d2 - d2 )l/2 ■ 

' TPAWQ'' 


( 20 ) 
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5. Initial conditions 

Since in both situations (low wavelength perturbations that os¬ 
cillate and long wavelength perturbations that grow) the initial 
state is 


p(gO = Po 
v(r,f) = 0, 



( 21 ) 

( 22 ) 


we set up the initial conditions as follows. 

First, A4 otal particles are distributed randomly within a 
unit cube and settled using non-self-gravitating SPH and peri¬ 
odic boundary conditions. This reduces the Poissonian density 
fluctuations, and produces an approximately uniform, but non¬ 
crystalline, density distribution (sometimes described as glass¬ 
like). The mean smoothing length of a particle is given by 


h = 


I 3 

\32;rAf„, 



1.1427 


(23) 


where we have substituted - 50, and the length of the 
edges of the cube is equal to unity. 

Second, a sinusoidal density perturbation is imposed by ad¬ 
justing the unperturbed .r-coordinate, x., of each particle i to a 
perturbed value, x', satisfying 


AA . 
— sin 


In 



Xj . 


(24) 


^/^JEANS 



This equation must be solved numerically for x' = x'(X;). We 
use a perturbation with fractional amplitude A = 0.1. 

Since periodic boundary conditions are being invoked, we 
can only apply perturbations which fit an integer number of 
wavelengths into the side of the unit cube, i.e. 

T = , (25) 


where = 1, 2, 3, etc.. 

A convenient measure of the resolution is the ratio of the 
mean diameter of an SPH particle, d - 4h, to the wavelength 
of the perturbation, A - i.e. 



Ah I ^ '^NEIB 

n. 4rt = n 


7tN„ 


(26) 


^/^JEANS 

Fig. 2. Characteristic timescales for the evolution of plane- 
wave perturbations, as a function of wavelength. The ordinate 
is the wavelength in units of the Jeans length, and the abscissa 
is the timescale in units of (Gp„)"^^^. For perturbations which 
oscillate (i.e. those with A < Tj^ans) the oscillation period esti¬ 
mated from the SPH simulations is represented by an open cir¬ 
cle. For perturbations which collapse (i.e. those with A > /Ijeans) 
the time for the peak density in the SPH simulations to increase 
by a factor 1.54 (see text) is represented by an open star. For 
reference, the analytic timescales are given by solid curves, (a) 
The very well resolved case, 7? = 0.5. (b) The marginally re¬ 
solved case, - 1.0. (c) The under-resolved case, "R - 2.0. 
Note that in all cases, even the under-resolved one, wavelengths 
which should oscillate do oscillate (i.e. no Jeans-stable pertur¬ 
bations artificially collapse). In (c), the filled circles and filled 


Thus a small value of R corresponds to good resolution, and 
the Jeans condition (Eqn.|3l can then be rewritten in the form 

R < 1. (27) 


where A = Tjeans ■ The number of SPH particles in one Jeans 
mass is then 








(28) 


The Jeans wavelength can be varied arbitrarily by changing 
the isothermal sound speed, a. 
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6. Test results 

In Figure ^ we compare the sinusoidal perturbation with the 
smoothing kernel for three values of the resolution: (a) “R - 0.5, 
i.e. very well resolved with = 400 SPH particles in one 

Jeans mass; (b) R - 1.0, marginally resolved, with - 50 
SPH particles in one Jeans mass; and (c) R - 2.0, under¬ 
resolved, with just A/fjEm/S - 6 SPH particles in one Jeans 
mass. We see why the resolution R - 1.0 is critical, in the 
sense that the smoothing kernel is closely matched to the per¬ 
turbation (see Fig.[0, and there are just enough SPH particles 
to resolve a three-dimensional object. 

In Figure|2]we plot the results of SPH simulations of plane- 
wave perturbations having different values of Each 

panel of Fig. |2l corresponds to a different resolution, R, viz. 
(a) R - 0.5, (b) R - 1.0, and (c) R = 2.0. If the perturbation 
oscillates, we plot with an open circle the oscillation period. 
If the perturbation grows, we plot as an open star the time re¬ 
quired for the amplitude to increase by a factor cosh(l) = 1.54. 
The analytic predictions for these times (Eans.[T^and l20> are 
shown as solid lines. Erom these plots, we can draw the follow¬ 
ing two important conclusions. 

(i) There is no artificial fragmentation. Even when the res¬ 

olution is poor (i.e. R - 2.0, EigEIc)), perturbations which 
should oscillate (d < do oscillate, and perturbations 

which should grow (d > do grow. Eurthermore we 

should be mindful that R - 2.0 corresponds to a strong vio¬ 
lation of the Jeans Condition, with only ~ 6 particles in a Jeans 
mass. It would never be tolerated in a simulation of collapse 
and fragmentation. 

(ii) With R ~ \ (good to marginal resolution), the char¬ 

acteristic timescales are reproduced well, except for the cases 
where d ~ . Marginally Jeans stable perturbations oscil¬ 

late a little too rapidly, and marginally Jeans unstable perturba¬ 
tions grow rather more slowly than they should. 

(iii) The effect of large R (poor resolution) is to shift the 
asymptote separating stable wavelengths from unstable ones, 
from d = djE,^^s to a slightly longer wavelength, i.e. to stabilise 
marginally Jeans unstable wavelengths, as if the temperature 
had been increased slightly. This is the reason for the gap on 
Fig. EJc) between d/djE,,Nj. = 1.0 and d/djE,,^;. = 1.4. A pertur¬ 
bation with (say) d = 1.2djE^Eis stabilised and initially oscil¬ 
lates, but at the same time, because of the poor resolution, its 
energy is transferred to longer wavelengths modes which then 
become unstable with a shorter timescale than the initial pertur¬ 
bation, and therefore it is not possible to evaluate an oscillation 
period. Thus poor resolution simply suppresses the collapse of 
marginally unstable modes. 

7. Correcting for the self-gravity of an individual 
SPH particle 

In standard self-gravitating SPH, the mutual gravitational force 
between two different SPH particles is included in the equation 
of motion, but the self-gravity of an individual particle is not. 
When the Jeans mass is not well resolved, we can improve the 
performance of the code by correcting for the fact that part of 
the pressure of an SPH particle must be used to support the 


particle against its own self gravity , rather than pushing on 
other particles. To formulate this correction, consider particle i 
in isolation. If its mass is m. and its sound speed is a., then 

3 J PdV = 3OT,af, (29) 

where the integral is over the volume of the SPH particle. We 
use a. because in general each SPH particle has its own sound 
speed which differs from that of other particles. The Virial 
Theorem tells us that, if the particle is to be in hydrostatic 
equilibrium, this integral must equal the magnitude of its self- 
gravitational potential energy, which is 

GirP'W 

|Q| = —(30) 
h. 


Here h. is its smoothing length, and 

r>s-2 ps'=s 

W = j ■I W(s')47r s'^ ds' ■ W(s)47rsds (31) 
Js=0 Js'=0 

is an integral which can be worked out given the dimensionless 
kernel function VT(i). It follows that the effective sound speed 
squared, should be reduced below the true sound speed 
squared, af, viz 


Gm.W 

3h. 


Eor the standard M4 kernel which we use, 

if 1 - I if 0 < s < 1; 

lV(i) = - J i(2 - i)3 if 1 < i < 2 ; 

^ [ 0 if s >2. 


(32) 


(33) 


(Monaghan & Lattanzio, 1985), W - 0.505 . 

We have repeated the simulations with R = 2.0, invoking 
this correction factor, and the results are presented as filled cir¬ 
cles and filled stars on Pig. Etc). As expected, the collapsing 
wavelengths collapse more rapidly (although still more slowly 
than the analytic solution). 


8. Conclusions 

The conclusions are very simple. SPH using the standard M4 
kernel and kernel-softened gravity (i.e. the standard options) 
only captures fragmentation which is (a) genuine, and (b) re¬ 
solved. It does not suffer from artificial fragmentation. 
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